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aijf  bijr 


NOMENCLATURE 


(See  Table  1  for  definitions  relevant  to  data 
input  quantities) 


*i 


Bi 


<*g/Pgo 

ac/ p CO 
■LcA.PgXg/Cg 
0 

Bi/Le 

ac^c/^c 


;  i  =  1-9,  13  (gas) 

:  i  as  10-12,  14  (solid) 


:  i  =  1-9 
:  i  =  10-12 
:  i  -  13 
:  i  *  14 


cij  ’  dij 
fli  (*j> 


Defined  in  Eq*  B4-B7 

AiBi  (1  +  tj)2  :  solid 

AiBi  (1  -  4 j) 2  •*  ?as 


,Ai  (1  +  tj)  *  solid 

f2i(tj)  (l  -  jj)  :  gas 


gi(us) 


jms  "  Bi 
\ms  +  Bi 


:  solid 
:  gas 


hi(tfj)  Volumetric  source  of  quantity  i  at  node  j 
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1.  INTRODUCTION 


The  governing  equations  of  a  unifird  theory  of  solid  propellant  Ignition 
were  developed  (see  Footnote  1)  and  are  summarized  in  Appendix  A.  Because 
of  the  occurrence  of  the  nonlinear  source  terms  V  and  S  and  of  the  surface 
regression  ms  there  is  little  chance  for  the  existence  of  an  exact  analytic 
solution  of  the  general  set  of  equations.  Simple  sub-theories  have  been 
successfully  handled  by  approximate  methods  such  as  matched  asymptotic  ex- 
pansions^  and  local  similarity3;  it  is  not  unlikely  that  other,  as  yet  un¬ 
treated,  sub-theories  may  eventually  yield  to  comparable  techniques. 

There  are  two  general  approaches  to  the  acquisition  of  solutions  to  the 
unsolved  sub-theories:  (1)  development  of  special  methods,  either  numerical 
or  approximate  analytic,  for  handling  each  sub-theory  in  an  optimum  manner 
and  (2)  development  of  a  generalized  numerical  program.  The  latter  course 
is  followed  here.  Even  though  this  approach  leads  to  a  much  more  compli¬ 
cated  computer  program,  it  has  the  advantage  of  greater  flexibility. 


2.  CHOICE  OF  METHOD  OF  SOLUTION 


In  addition  to  the  approximate  analytic  techniques  already  mentioned, 
there  are  a  great  number  of  numerical  methods  for  obtaining  solutions  to 
partial  differential  equations.  These  methods  fall  generally  into  two 
classes:  (1)  direct  replacement  of  the  derivatives  by  appropriate  finite 
difference  analogs  and  (2)  use  of  an  indirect  approach,  typified  by  the 
methods  of  weighted  residuals, which  often  leads  to  a  finite  difference 

1  Naval  Weapons  Center.  A  Unified  Theory  of  Solid  Propellant  Ignition ; 
Part  l  -  Development  of  Mathematical  Modelt  H.  H.  Bradley,  Jr.,  China  Lake, 
CA,  NWC  August  1974  (NWC  TP  5618,  Part  1). 

2  Linan,  A.,  and  F.  A.  Williams,  "Radiant  Ignition  of  a  Reactive  Solid 
with  In-Depth  Absorption,"  COMBUS  AND  FLAME,  Vol.  18  (1972),  pp.  85-97. 

3  Waldman,  Cye  H.  "Theory  of  the  Heterogeneous  Ignition,"  COMB  SCI 
TECH,  Vol.  2  (1970),  pp.  81-93. 

*  Finlayson,  B.  A.,  and  L.  E.  Scriven.  "The  Method  of  Weighted  Residuals 
—A  Review,"  APPL  MECH  REV,  Vol.  19,  No.  1  (September  1966),  pp.  735-48. 

5  Crandall,  S.  H.  Engineering  Analysis ,  A  Survey  of  Numerical  Pro¬ 
cedures.  New  York,  McGraw-Hill  Book  Co.,  1956. 
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formulation.  The  latter  results  In  a  smoothing  effect  which  may  tend  to 
obscure  details  of  solutions  varying  rapidly  with  position  unless  many 
weighting  functions  are  used.  In  addition,  the  algebraic  manipulations  be¬ 
come  somewhat  tedious  for  large  systems.  In  applying  the  direct  finite 
difference  methods,  care  is  essential  to  avoid  numerical  instabilities 
which  often  obscure  the  true  solution. 

Three  possibilities  exist  in  the  direct  application  of  finite  difference 
techniques  to  the  set  of  parabolic  equations  of  the  unified  ignition  theory 
depending  upon  whether  one  or  both  independent  variables  (time  and  space) 
are  subjected  to  finite  difference  operations. 6  Differencing  in  the  time 
direction  (often  referred  to  as  the  method  of  lines)  produces  a  large  system 
of  simultaneous  ordinary  differential  equations  which  are  stiff  by  nature. 

Some  special  methods,  such  as  Gear's^  are  available  for  solution  of  such 
systems;  good  stability  and  accuracy  with  reasonable  computing  times  are 
bought  at  the  expense  of  very  large  computer  requirements.  The  final  method, 
differencing  along  both  the  time  and  space  axes,  can  be  carried  out  either 
explicitly  or  implicitly.  In  most  explicit  methods,  computing  parameters 
are  severely  limited  owing  to  stability  considerations,  resulting  in  excess¬ 
ive  computing  times.  In  general,  implicit  methods  are  more  stable  but  re¬ 
quire  the  solution  of  simultaneous  nonlinear  algebraic  equations;  nevertheless, 
the  special  form  of  the  system  matrix  after  linearization  through  iterative  or 
predictor-corrector  methods  makes  the  implicit  approach  more  attractive. 

Two  versions  of  the  implicit  method  will  be  described,  differing  principally 
in  the  treatment  of  the  nonlinear  terms.  The  computer  code  is  implemented 
for  one  method  and  may  easily  be  extended  to  the  ether  version. 


3.  DESCRIPTION  OF  COMPUTING  ALGORITHMS 


3.1  FIELD  EQUATIONS 


The  equations  to  be  solved  are  displayed  in  Appendix  A.  For  convenience 
in  describing  the  computing  algorithm,  the  field  equations  (Eq.  A1  -  A2)  may 
be  written  compactly  as: 


3v 

3t 


d2v 

-7  + 
1  9<J> 


9v 

9<J> 


*  h(v) 


(1) 


6  Crank,  J.,  and  P.  Nicholson.  "A  Practical  Method  for  Numerical  Evalu¬ 
ation  of  Solutions  of  Partial  Differential  Equations  of  the  Heat-Conduction 
Type,"  CAMBRIDGE  PHIL  SOC,  PROC,  MATH  PHYS  SCI,  Vol.  43  (1947),  pp.  50-67. 

7  Gear,  C.  W.  Numerical  Initial  Value  Problems  in  Ordinary  Differential 
Equations .  Englewood  Cliffs,  N.  J.  Prentice— Hall,  1971.  253  pp. 
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In  Ea.  (1),  v  represents  the  vector  dependent  variable  with  mass  fraction 
and  enthalpy  as  scalar  components;  fj  and  f2  are  functions  o  pos  ^ 

(nodal  index)  only  which  arise  from  the  co°fdinft«  J t!®"erf^  and  h ±8*1 
is  a  nonlinear  function  of  v  evaluated  at  Jsing 

nonlinear  function  of  v  (hence  of  4)  evaluated  in  the  bulk We.  u  g 
central  differences  to  approximate  the  first  and  second  derivatives 

obtain: 


/  k+1 

\Uj 


»5)/At  ■  wJK1  *  a'9’“j] 


+  f. 


.wsK1 


(l-0)u°  g 


(».)  *  h(“j) 


(2) 


(3) 


(4) 


where  the  difference  operators  are  defined  as: 

£  (VI  ■ 

\wj  -  Ki  •  wj-i  Ym 

Equation  (2)  represents  a  8*""*1  “^“^tacripMTiS 

solution  from  tima  t  through  th.  int«v.l  it.  Th.^  P  ^  }  correB_ 

denote  approximations  to  the  tru  *lvll  The  terns  g  and  S  represent 

ponding  to  the  times  t  and  t  ,  P  interval  At.  The  parameter 

appropriately  chosen  values  ° ritta^lB  fully  implicit  (backward  difference 
6  determines  whether  the  algorithm  values  of  9  may  be  chosen, 

algorithm,  9=1)  or  explicit  (  )•  .  ,  scheme,  being  the  most  common. 

0-1/2,  corresponding  to  the  Crank-Nicholson  scheme,  n  8  to 

The  explicit  scheme,  with  *  and  * b '  -aluated^at  ^curren^  ^  ^  ^i; 

linear  algebraic  equations  w|)*ch  restrict  its  use.  All  other 

however,  stability  considers  o  Q  >  jn)  generate  systems  of  simul- 

schemes  (stability  is  possi  e  on.^.  ,  “  Qf  uk+l  at  three  nodes  per 

taneous,  algebraic  discussed  later).  The  nonlinear- 

equation  (except  at  the  boundaries,  be  han(jied  by  any  of  a 

ities  of  the  system  in  the  present  baronies  o f  which  are  described  in  the 
number  of  linearization  techniques,  examples  - 

following  sections. 

3.1.1  Iterative 
EH.  parameter  .. 

STS.’SSSS-  £££  at  the  current  and  advanced  times.  It  seems 
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logical  to  apply  a  weighting  function  to  the  nonlinear  functions  g  and  h. 
If  the  same  weighting  function  la  used,  then 


9(us)  *  ?[e“s  * 

=  h[eu*  +  a-B)u°  j 


(5) 


(6) 


where  the  superscript  *  denotes  the  previous  approximation  to  the  unknown 
function  at  t  +  At.  It  may  be  obtained  by  extrapolation  from  earlier  values 
of  the  solution  or  by  simply  using  the  value  of  the  function  at  t  as  a 
first  approximation.  The  resulting  iterative  process  (Picard's)  is  linearly 
convergent.  A  more  sophisticated,  quadratically  convergent  iteration 
(Newton's)  involves  a  series  expansion  of  the  functions  u£  and  with  the 
linear  terms  being  retained.**  In  the  current  application^  it  is  believed 
that  the  number  of  evaluations  of  complicated  functions  required  by  this 
latter  method  would  increase  computing  times  more  than  is  warranted. 


3.1.2  Predictor-Corrector 


The  algorithm  actually  utilized  in  the  computer  coding  consists  essen¬ 
tially  of  a  Crank-Nicolson  step  of  At  with  nonlinear  functions  of  dependent 
variables  evaluated  at  At/2.  These  latter  values  are  first  estimated  by  an 
implicit  predictor  step  (8 *1)  of  A t/2.^  We  thus  have  for  the  predictor: 


k  =  0,  0  »  2,  At  -*■  At/2 
and  for  the  Crank-Nlcholson  corrector 
*-2,0-’  2/2,  At  -*•  At 


where 

g(us)  =  g(uks) 

h( Uj)  -  h(uk) 


8 

Ames,  William  F.  numerical  Methods  for  Partial  Differential  Equations. 

London,  Thomas  Nelson  and  Sons,  Ltd.,  1969. 

9 

Douglas,  J.,  and  B.  F.  Jones.  "On  Predictor-Corrector  Methods  for  Non- 
Linear  Parabolic  Differential  Equations,  SOC  1ND  APPL  MATH,  J.  Vol.  11,  No.  1 
(March  1963),  pp.  195-204. 
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When  the  above  relationships  are  substituted  Into  Eq 
becomes 


(2),  the  predictor 


1 

u  ,  - 


it/2 


'  2i 


and  the  corrector. 


(7) 


0\ 

u4 


/2  + 


u°\gr ./  u1  \/2  +  h 
jril  si 


(8) 


3.2  INTERFACE  EQUATIONS 

The  matching  conditions  at  the  pha.elnterf.ee  (Eq.  ««-“> 
derivatives  of  the  unknowns  (excep  s°  Using  central  differences  at 

interface  a.  well  as  value.  “""”°™a.eo"s  and  solid  nodal  points.  If 

rra^dThic  r°fiueid  «. 

SSTfiSd*1^^1  -  'hctUlous  values  me, 

solved^for  ^  ^  ^  “ 

shown  in  Appendix  B . 

Special  treatment  is  require  for 

derivative  of  equation  by  which  fictitious  nodal 

ZtSZ  be* eliminated  as  In  the  KME- 

In  the  application  of  the  Predator-corrector  ®1®J'1£|!"nMIlchoi,on  form 
equations,  a  backward  ^  S'^'mtU.  which  can  be  in- 

~  tL  precautlon  ”s  taun 
until  further  investigations  could  be  completed. 


10 


-  Douglas,  J.  "A  Survey  of  Numerical  Methods  for  Parabolic  Differential 
Equations,"  in  ADVANCE  COMPUTER ,  New  York,  Academic  Press,  1961.  Vol  , 
pp.  1-54. 
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4.  DESCRIPTION  OF  COMPUTER  PROGRAM 


The  computer  program  was  designed  to  take  advantage  of  the  capabilities 
of  the  Univac  1108  Exec  8  system  at  the  Naval  Weapons  Center.  Coding  is  in 
FORTRAN  V  language,  use  being  made  of  PARAMETER  and  NAMELIST  statements  and 
internal  subprograms. 

The  input  to  the  computer  program  consists  of  the  assignments  for  116 
quantities  representing  physical  parameters,  program  controls,  and  output 
specifications.  To  simplify  program  use,  standard,  or  base  values  are 
assigned  to  80  of  these  quantities  (see  Table  1) ;  other  than  base  values 
are  introduced  by  way  of  a  NAMELIST  statement.  The  remaining  36  quantities 
are  specified  through  conventional  input  methods. 

The  solution  of  the  equations  includes  values  of  the  fourteen  dependent 
variables  at  each  time  interval  and  each  space  node  of  the  computing  mesh. 
Since  a  minimum  of  100  nodes  has  been  selected  to  assure  accuracy,  the  total 
output  of  a  run  of  100  to  300  time  steps  would  represent  an  overwhelming 
volume  of  printout.  As  a  compromise,  only  a  fraction  of  the  computed  output 
is  printed  during  the  initial  program  execution;  selection  of  the  printed 
portion  is  controlled  by  input  parameters.  To  avoid  loss  of  the  unprinted 
solution,  the  entire  solution  is  stored  on  tape,  through  intermediate  use 
of  high  speed  mass  storage,  for  later  processing  and  analysis,  which  may  in¬ 
clude  plotting  or  additional  printing  of  results. 


4.1  FLOW  CHART 

A  simplified  flow  chart  of  the  program  is  shown  in  Figure  1.  A  complete 
listing  of  the  program  is  available.  Address  Commander,  Naval  Weapons  Center, 
Code  608,  China  Lake,  CA  93555.  The  contents  of  each  block  is  briefly  pro¬ 
vided  in  the  following  descriptions. 

PROGRAM  INITIALIZATION  -  Contains  specification  statements  such  as  array 
dimensions,  formats,  data  statements  (including  specification  of  the  values 
of  the  base  set  of  80  parameters),  and  several  program  constants.  The  out¬ 
put  mass  storage  file  defined  into  which  the  complete  solution  is  to  be 
stored  (ISF=10) .  This  unit  number  may  be  changed  to  conform  to  local  com¬ 
puter  requirements.  In  the  current  coding,  only  one  such  file  may  be  gene¬ 
rated  per  program  execution. 

RUN  INITIALIZATION  ~  (1)  Saves  the  standard  parameter  set  for  later  res¬ 
toration.  This  provision  is  included  to  allow  several  run.;  to  be  stacked 
and  would  require  defii  ition  of  an  equivalent  number  of  storage  files.  (2) 
Uses  NAMELIST  input  to  alter  selected  base  parameters;  uses  standard  input 
(READ)  statements  for  specifying  run  identification  and  output  print  con¬ 
trols. 


8 


NWC  TP  5618,  Part  2 


TABLE  1,  NAMELIST  Input  Parameters 
and  Standard  Values. 


Computer 

Algebraic 

Standard 

Definition 

designation 

designation 

value 

FREQ  1 

V1 

0. 

Pre-exponential  factor,  reaction  10a 

FREQ  5 

v5 

0. 

Pre-exponential  factor,  reaction  5b 

FREQ  6 

v6 

0. 

Pre-exponential  factor,  reaction  6C 

FREQ  7 

v7 

0. 

Pre-exponential  factor,  reaction  7d 

FREQ  8 

v8 

0. 

Pre-exponential  factor,  reaction  8e 

FREQ  10 

v10 

0. 

Pre-exponential  factor,  solid  fuel 

vaporization 

FREQ  11 

V11 

0. 

Pre-exponential  factor,  solid  oxidizer 

vaporization 

FREQ  12 

v12 

0. 

Pre-exponential  factor,  solid  product 

vaporization 

El 

El 

0. 

E5 

E5 

0. 

E6 

E6 

0. 

Activation  energy.  Subscripts 

E7 

E7 

correspond  to  those  for  pre- 

U. 

exponential  factors. 

E8 

E8 

0. 

E10 

ElO 

0. 

Ell 

E11 

0. 

E12 

e12 

0. 

Q1 

Qi 

0. 

Q5 

Q5 

0. 

Q6 

Q6 

0. 

Q7 

Q7 

0. 

Heat  of  reaction  (positive  for 

Q8 

Q8 

exothermic) .  Subscripts  corres- 

u. 

pond  to  those  for  pre-exponential 

Q10 

Qio 

0. 

factors. 

QH 

Qn 

0. 

Q12 

Ql2 

0. 
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TABLE  1.  (Contd.) 


Computer 

designation 

Algebraic 

designation 

Standard 

value 

VF7 

vf7 

1. 

VO  7 

vo7 

1. 

VF8 

VF8 

1. 

V08 

V08 

1. 

SOS 

S5 

1. 

S06 

S6 

1. 

STOV2 

b2 

1. 

ST0V4 

b4 

1. 

STOV  11 

bll 

1. 

ST0S2 

b2 

1. 

STOS4 

b4 

1. 

MU1 

1. 

MU2 

v2 

1. 

MU4 

”4 

1. 

MU  5 

1. 

MU  6 

*6 

1. 

MU  7 

V7 

1. 

MU8 

1. 

MU9 

w9 

1. 

MU10 

v10 

1. 

MU11 

un 

1. 

MU12 

1. 

COND  (1) 

X 

c 

1. 

COND  (2) 

*8 

1. 

SPHT  (1) 

Cc 

1. 

SPHT  (2) 

C8 

1. 

Definition 

Fuel  reaction  order,  reaction  7 
Oxidizer  reaction  order,  reaction  7 
Fuel  reaction  order,  reaction  8 
Oxidizer  reaction  order,  reaction  8 
Oxidizer  reaction  order,  reaction  5 
Oxidizer  reaction  order,  reaction  6 
Oxidizer  stoichiometry,  reaction  7 
Oxidizer  stoichiometry,  reaction  8 
Oxidizer  stoichiometry,  reaction  10 
Oxidizer  stoichiometry,  reaction  5 
Oxidizer  stoichiometry,  reaction  6 
Mol.  wt.  of  vaporized  solid  fuel^ 

Mol.  wt.  of  vaporized  solid  oxidizer 
Mol.  wt.  of  original  gaseous  oxidizer 
Mol.  wt.  of  product  in  reaction  5 
Mol.  wt.  of  product  in  reaction  6 
Mol.  wt.  of  product  in  reaction  7 
Mol.  wt.  of  product  in  reaction  8 
Mol.  wt.  of  inert  gaseous  diluent 

Mol.  wt.  of  solid  fuel 

Mol.  wt.  of  solid  oxidizer 

Mol.  wt.  of  solid  product 

Solid  thermal  conductivity 
(cal/cm  sec  deg  K) 

Gas  thermal  conductivity 

Solid  specific  heat 
(cal/gm  deg  K) 

Gas  specific  heat 
(cal/gm  deg  K) 
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TABLE  1.  (Contd.) 


Computer 

designation 


Algebraic 

designation 


Standard 

value 


RH01 

pco 

1. 

T.E 

Le 

1. 

BE 

8 

1. 

TLIM(l) 

T 

Co 

1. 

TLIM(2) 

T 

go 

1. 

PRESS 

p 

82.  C 

YFZ1 

Yfc 

1. 

Y0Z2 

YOg 

1. 

QDOT 

qr 

1. 

TMELIM 

10. 

TSLIM 

3. 

FCT 

2. 

NREFIR 

0. 

TCO 

11. 

CUTOFF 

F 

SUB 

F 

Definition 


Solid  density  (gm/cm3) 

Lewis  number  of  gas  •  p cV/\ 

Solid  extinction  coefficient  cm_l 
Initial  solid  temperature  °K 
Initial  gas  temperature  °K 
Pressure  (atmospheres) 

Initial  gas  fuel  mass  fraction® 
Initial  gas  oxidizer  mass  fraction 
Radiant  flux  (cal/cm2  sec) 

Ignition  time  limit  (sec) 

Ignition  temperature  limit  (deg  K) 

One  greater  than  number  of  heating 
periods  allowed  for  ignition  after 
cutoff 

Number  of  refires  allowed  to  achieve 
ignition  in  cutoff  case 

Time  to  first  cutoff 

Cutoff  flag:  F  -  no  cutoff 
T  ■  cutoff 

Subsurface  absorption  flag: 

F  ■  no  SSA 
T  -  SSA 


REG 

F 

Surface  regression  flag: 

F  -  no  regression 

T  ”  regression 

ALP(l) 

ac 

1. 

X-scale  transformation  for  solid 

ALP (2) 

°« 

1. 

X-scale  transformation  for  gas 

N(l) 

Nc 

NNN-1 

Number  of  strips  in  solid 

N(2) 

N8 

NNN-1 

Number  of  strips  in  gas 

'-iA 
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TABLE  1.  (Contd.) 


Computer 

designation 

Algebraic 

designation 

Standard 

value 

Definition 

DULIM 

.01 

DULIM  T/T  in  search  for  time  step 

DLLIM 

.005 

DLLIM  T/T  in  search  for  time  step 

SPR05 

0 

Flag  for  Including  reaction  5* 

SPR06 

0 

Flag  for  including  reaction  6* 

SPR010 

0 

Flag  for  including  solid  fuel 

pyrolysis* 

SPR011 

0 

Flag  for  including  solid  oxidizer 

pyrolysis* 

SPR012 

0 

Flag  for  including  solid  product 

pyrolysis* 

NHLIM 

10 

Number  of  time  step  halvings  allowed 

in  search  for  time  step 

SURF 

F 

Flag  for  all  surface  processes 

r  -  ignore,  T  -  include 

a  Solid  Reaction:  [F]  +  b^o]  -  b12[P] 
b  Surface  Reaction:  [F]  +  bjlo]  *  b5[P] 
c  Surface  Reaction  [F]  +  b^[oJ  -  bg [P ] 
d  Gas  Reaction:  [F]  +  b2[o]  -  b7 [P] 
e  Gas  Reaction:  [F]  +  b^fo]  -  bg[P] 

f  Set  all  molecular  weights  equal  in  this  version  of  the  computer  program 

g  Set  Yfc  to  any  negative  value  to  obtain  reverse  of  reaction  10.  Then 
the  program  will  set  Ypc  -  1.  and  Yfc  -  YQ  -  0.  Also  vj  must  be  set 
negative,  and  Qj  becomes  negative  for  exothermic  reaction. 


' <• juals  0,  do  not  include  process;  equals  1,  include  process. 
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(3)  Prints  and  (optionally)  stores  the  run  identify lcat: ion  and  parameters. 

(4)  Calls  subroutine  START  to  compute  run  constants  from  input  alters 
and  to  estimate  initial  time  step  for  integration.  (5)  Calls  wbrout ine 
TSUB  to  convert  enthalpies  to  temperatures.  Prints  and  (optionally)  stores 
initial  conditions.  (6)  Calls  subroutine  ANFORM  to  evaluate  solution  at  end 

of  first  time  step  by  analytical  formula  from  which  n°nli™^  PRCS  is  called 
effects  of  subsurface  absorption  are  omitted.  The  subroutine  ERCS  call 
to  evaluate  the  error  functions  occurring  in  the  analytic  solution.  It  al 
prints  and  (optionally)  stores  the  analytic  step. 


ADVANCE  SOLUTION  FROM  t  TO  t  +  At-  This  segment  of  coding  contains 
(through  subroutines  ABCCOE,  DCOE,  TRIDAG,  SPEC,  S,  VOL,  AND  S1314)  the 
implicit/Crank-Nicholson  predictor-corrector  algorithm.  Program  flow  is 
directed  in  accordance  with  input  options. 


TEST  ACCEPTABILITY  OF  SOLUTION  -  Subroutine  CRIT  contains  the  coding 
which  ascertains  the  acceptability  of  the  solution  after  each  time  step. 
Criteria  examined  are:  (1)  non-negativity  of  concentration  and  (2)  bracket¬ 
ing  of  the  time  rate  of  change  of  temperature  between  upper  and  lower  abso¬ 
lute  limits.  If  a  criterion  is  violated  the  time  step  is  halted  or  doubled 
as  appropriate  and  the  step  repeated.  This  step  size  adjustment  is  not  in¬ 
voked  until  after  the  first  ten  steps  and  is  then  limited  to  a  number  of 
halvings  as  specified  by  an  input  control  parameter.  When  this  number  is 
equalled  the  run  is  terminated.  In  addition,  CRIT  controls  the  go/no-go 
computations  as  follows: 

(a)  Constant  heat  flux  is  applied  until  a  specified  time,  when  the  flux 
is  reduced  to  zero.  The  solution  at  that  time  is  saved. 

(b)  The  solution  is  continued  with  zero  flux  for  a  specified  time  inter 
val  after  shutoff. 


(c)  If  the  temperatures  at  the  surface  and  at  the  first  ten  gas  nodes 
are  decreasing  after  the  above  specified  time  interval,  then  the  saved  cu  - 
off  values  are  restored,  and  heating  is  continued  for  one  more  time  step  be¬ 
fore  another  shutoff  is  invoked.  Steps  (b)  and  (c)  are  then  repeated.  The 
search  for  ignition,  as  evidenced  by  an  increase  in  temperature  at  the  sur¬ 
face  or  at  one  or  more  of  the  first  ten  gas  nodes,  is  continued  for  a  speci¬ 
fied  number  of  trials  before  the  entire  run  is  terminated. 


OUTPUT  SOLUTION  -  Prints  solution  after  a  specified  number  of  time  steps 
and  space  nodes  (not  necessarily  equal  increments)  P.nd  stores  (optionally) 
the  entire  solution  after  every  time  step. 


END  OF  RUN  -  When  specified  heating  time  or  surface  temperature  is 
reached,  the  base  parameter  set  is  restored  and  a  new  run  may  be  started. 
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4.2  PROGRAM  INPUT 

(a)  First  two  cards  (Columns  2-80).  Any  message,  identification,  or 
special  notes  concerning  the  run.  The  information  contained  is  stored  on 
the  mass  storage  file  when  the  appropriate  option  is  used,  and  would  be 
valuable  in  the  process  of  file  searching  during  data  retrieval. 

(b)  As  mentioned  in  earlier  sections,  the  program  in  its  current  form 
requires  the  specification  of  116  quantities  for  its  execution.  To  facili¬ 
tate  the  preparation  of  input,  standard  values  are  assigned  (through  a  DATA 
statement)  to  80  quantities,  which  are  listed  in  Table  1.  By  the  use  of  the 
NAMELIST  Provision,  any  combination  of  values  of  the  80  quantities  may  be 
modified  by  simply  specifying  the  new  values.  Details  of  using  NAMELIST 
are  available  at  computer  sites  which  utilizes  a  FORTRAN  V  compiler.  If 
NAMELIST  is  not  available,  the  input  statements  must  be  modified  by  sub¬ 
stituting  conventional  coding  with  formatted  READ  statements. 

(c)  Final  card  (format  1411,  IX,  1513,  611,  F6.0)  used  with  a  conven¬ 
tional  READ  statement  to  input  the  following  controls: 

ICALC  (Columns  1-14)  -  a  fourteen-element,  single-digit  integer  vector 
which  determines  which  of  the  dependent  variables  are  to  be  computed.  When 
ICALC (I)  is  set  equal  to  zero,  calculations  for  dependent  variable  I  are  by¬ 
passed,  otherwise  the  calculation  is  made.  ICALC (14)  is  set  to  one  by  the 
program,  therefore,  solid  phase  temperatures  are  always  calculated. 

NC  (Columns  16-36)  -  a  seven-element,  array  specifying  the  solid  phase 
nodes  to  be  printed  out.  The  nodes  are  counted  from  the  surface  as  number 
one  and  increase  leftward  into  the  solid  phase.  No  decimal  point  should  be 
punched.  Right  justify  the  integers  in  the  appropriate  fields.  No  integer 
exceeding  the  total  number  of  solid  nodes  should  be  punched  (101  as  currently 
programmed) . 

NG  (Columns  37-57)  -  Similar  to  NC  except  for  gas  phase  nodes.  Counting 
is  from  the  surface  as  the  number  one  gas  node  (hence,  duplicating  the  number 
one  solid  node)  and  proceeding  rightward.  The  same  precautions  apply. 

NPINT  (Columns  58-60)  -  Frequency  of  printout  (timewise).  NPINT  ■  0  or 
1  prints  every  step;  NPINT  ■  10  every  tenth  step,  etc. 

NCDENS  (Column  61)  -  Control  specifying  mode  of  computing  gas  phase 
density. 

OUTST  (Column  62)  -  Control  Integer  to  cause  solution  to  be  stored  on 
unit  ISF.  (ISF  *  10  in  coding.)  Must  not  be  zero  or  blank  if  storage  is 
desired. 
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OUTPL  (Column  63)  -  Control  Integer  to  provide  on  line  printer  plotting. 
Invokes  subroutine  PLOTTT,  details  of  which  would  be  dependent  on  local 
facilities.  Since  PLOTTT  is  now  vacuous,  OUTPL  has  no  effect  and  m|ty  be 
left  blank. 

NER  (Columns  64-66)  -  Diagnostic  debug  print  controls.  Leave  blank. 

FMDT  (Columns  67-72,  Format  F6.0)  -  Factor  by  which  initial  estimated 
time  step  is  multiplied  in  attempt  to  achieve  starting  accuracy.  Now  used 
as  1.0.  Will  be  set  to  1.0  if  left  blank  or  assigned  a  value  less  than 
0.001 


4.3  PROGRAM  OUTPUT 

The  first  page  of  printed  output  contains  the  message  input  on  the  first 
two  cards,  followed  by  a  matrix  listing  of  the  80  NAMELIST  quantities  and 
the  program  controls.  On  the  succeeding  pages  the  solution  is  arranged  five 
time  steps  to  the  page.  A  sample  printout  is  shown  in  Appendix  C  and  is 
self  explanatory  except  for  the  abbreviations.  The  line  identifying  the 
time  also  gives  the  stimulating  energy  flux  and  average  surface  mass  fluxes 
for  the  time  it  *:erval  with  the  starting  value  TIME;  also  provided  is  the 
number  of  times  the  integration  interval  was  halved  to  satisfy  the  accep¬ 
tance  criterion  of  subroutine  CRIT.  The  first  line  following  the  time 
records  a  reduced  temperature  (,T-T0)/Tq  for  specified  nodes  of  the  solid 
and  gas  phases.  Succeeding  lines  contain  species  mass  fractions.  In  gen¬ 
eral,  solid  nodes  are  to  the  left  while  gaseous  nodes  are  to  the  right.  To 
conserve  space  three  of  the  gaseous  species,  as  noted,  are  displayed  on  the 
left  and  should  be  associated  with  the  gaseous  nodes.  Following  is  the  key 
to  the  abbreviations : 

F(C)  -  on  left,  solid  fuel  (I  -  10);  on  righ4  gaseous  pyrolyzed  fuel 

(I  -  1). 


0(C)  -  on  left,  solid  oxidizer  (I  -  11);  on  right,  gaseous  pyrolyzed 
oxidizer  (I  -  2) . 

P(C)  -  on  left,  product  of  solid  reaction  (I  -  12);  on  right,  gaseous 
pyrolyzed  product  of  solid  reaction  (I  ■  3) . 

P(S,  F+0(C) )  -  gaseous  product  of  surface  reaction  between  solid  fuel 
and  pyrolyzed  oxidizer  (I  ■  5) . 

P(S,  F+0(G))  -  gaseous  product  of  surface  reaction  between  solid  fuel 
and  initial  gaseous  oxidizer  (1*6). 

INERT  -  gaseous  inert  diluent  (I  ■  9)  . 
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0(G)  -  Initial  gaseous  oxidizer  (I  -  4) . 

P(G,  F+0(C) )  -  gaseous  product  of  gas  phase  reaction  between  pyrolyzed 
fuel  and  pyrolyzed  oxidizer  (I  -  7) . 

P(G,  F+0 (G) )  -  gaseous  product  of  gas  phase  reaction  between  pyrolyzed 
fuel  and  initial  gaseous  oxidizer  (I  -  8) . 
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Appendix  a 


GOVERNING  EQUATIONS* 


Field  Equations  (species  and  enthalpy) 


Gas  (i  -  1  -  9,13 ;  <j>  >  0) 


. 2  l 


I  0  V*  , 

r  -  Vi «-«  -rr  -  "-♦'(  vsi )  if  *  'j 


Solid  (i  =  10  -  12,  14s  Bi  =  0,  i  «  10  -  12;  <P  <  0) 


dZ.  -  32Z.  3Z 

-sr  ■  Vi  <***>  -^r  -  Ai  w  *  vi 


Interface 


Species  (i  =  1  -  9;  mi+9  =  0,  i  >  3) 


32  . 

m.,„  *  m  z.  -  fl.  xr  -  5. 
2+9  si  a.  3<f>  i 


Enthalpy 


3ZJ«  3zi3 

B14  3*  BIJ  3*  *  S1314 


Remote  Boundary  Conditions 


Gas  (i  =  1  -  9,  13) 


(l,t)  -  0 


Solid  (i  =  10  -  12,  14) 


(-l,t)  =  0 


Initial  Conditions  (1=1-  14) 


Z  .  =  Z  . 

1  10 


See  Ref.  1  for  Nomenclature. 


Preceding  page  blank 


(A-l) 


(A-2) 


(A-3) 


(A-4) 


(A-5) 


(A-6) 


(A-l) 
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Appendix  B 

REPRESENTATION  OP  GOVERNING  EQUATIONS  IN  FINITE  DIFFERENCE  FORM 


The  finite  difference  formulation  of  the  equations  of  Appendix  A  Is 
based  on  the  two  step  predictor-corrector  algorithm  described  In  Section 
3.  A  finite  difference  net  is  established  using  equal  spacing  in  the 
<Mpace.  The  number  of  strips  in  the  solid  and  gas  phases  is  Ba  and  Ng , 
respectively,  with  nodal  numbers  running  from  one  (solid  remote  boundary) 
to  i'k+l  at  the  interface  and  continuing  to  %+fyj+l  at  the  gas  phase  remote 
boundary.  The  surface  node  is  repeated  to  allow  for  surface  discon¬ 
tinuities  in  species  concentration.  The  spacing  of  the  nodes  isAAi-2/W 
and  A $2-2/)^.  1  c 


The  predictor  and  corrector  equations  follow  from  Eq.  (7-8)  and  take 
the  following  forms: 


+  f 


2i 


W[u; 


7+1 


(B-l) 


u2 

-L 


At 


o 
u  . 

-J.  = 


A  Tu? 


7+2 


2  2 

-  2u.  +  u 


+  u 


7"1  j+1 


2u 


*  Vi] 


/2 


+  M**)  [vi  -  uU  +  Vi  -  v]’i(“!)"  *  IB-2, 

When  unknowns  are  collected  on  the  left  and  known  quantities  on  the 
right,  the  following  set  of  equations  is  obtained: 


Where 


k  k+1  k  k+1 

a  u  +  b .  .u  . 
7-1  i 7  7 


k  k+1 

C  .  .U  . ,  , 

17  7+1 


'  J  +  fii(*;)At 

-  hMYt/2 


for  the  predictor  (k  -  0) . 


( B-3 ) 

(B-4) 

(B-5) 

(B-6) 

(B-7) 


Preceding  page  blank 
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The  equations  for  a,  b  and  c  are  unchanged  for  the  corrector  (k *1); 
the  expression  for  d  becomes: 


i  -  fn 


'  hi  (“?)“  *  [fU  (*))  -  Lt/!  *  [• 

’  pii  (*j)  *  f2i(*j)s(^)]vi Lt/2 


The  finite  difference  forms  of  the  interface  equations  are: 
Gas  species: 


N^fi1 

(B-8) 


m 


i+9 


k  k+1  /  k+1  k+1  \  ....  -k 

=  m  u .  -fl.lu.  ,  -  u .  ,  }/2A4>  -  S  . 

s  i,s  i\  i,s+l  i,s-l/  i 


(B-9) 


Elimination  of  us_j  between  Eq.  (B-9)  and(B-3)  evaluated  at  j=s  (surface) 
provides  the  appropriate  form  for  the  tridiagonal  matrix. 

Solid  species:  Use  the  backward  difference  formula  for  3u/3<J>  in 
Eq.  (A-2) 


(B-10) 


if’  (“i-2  - 

and  eliminate  us_2  by  use  of  Eq.  (B-3)  evaluated  at  j  =  s-1. 

Enthalpy: 

-  »U,s-l)^c  -  *12  (-5U  -  W.s-^9  *  <°-n> 


The  quantities  Ui4/S+j  and  u i3fS-l  represent  fictitious  nodes.  They  are 
eliminated  by  means  of  the  solid  enthalpy  field  equation  for  ujj  and  ujj 
centered  at  the  surface.  The  resulting  equation  contains  u23,s  an<*  u14,s 
related  by 

“13,S-C*‘114.S/*  +‘ilTc°  -%>)  IB-12> 

which  arises  from  continuity  of  temperature  at  the  surface. 

The  result  of  the  above  operations  is  a  set  of  tridiagonal  equations 
to  be  solved  for  each  of  the  twelve  species  and  enthalpy  for  each  half 
of  the  computing  algorithm.  There  are  %+l,  N  +l  and  Afc+AL+1  equations 
in  the  tridiagonal  sets  for  solid  species,  gadeous  species  and  enthalpy, 
respectively. 
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Appendix  C 

SAMPLE  PRINTOUT  OF  COMPUTER  CALCULATIONS 

The  following  two  pages  show  an  example  of  the  output  produced  by  the 
computer  program.  The  entire  run  is  not  reproduced. 

The  first  page  shows  a  matrix  listing  of  the  80  NAMELIST  variables  and 
their  associated  numerical  values.  The  particular  values  shown  for  FREQ  1, 
El,  Ql,  COND  (1),  SPHT  (1),  RHOl,  TLIM  (1),  YFZ1,  and  QDOT  were  chosen  to 
duplicate  a  previous  computation  for  solid  phase  ignition.  The  interpre¬ 
tation  of  the  final  line  on  the  first  page  is  given  in  Section  4.2  (Program 
Input)  paragraph  (c). 

The  second  page  shows  a  sample  output  for  the  first  five  time  steps. 

A  complete  explanation  is  given  in  Section  4.3  (Program  Output).  Since 
only  solid  phase  calculations  were  made  in  this  example,  the  results  per¬ 
taining  to  the  gas  phase  appear  as  zeros. 


Bradley,  H.  H. , 
Constant  Energy  Flux," 


Jr.  "Theory  of  Ignition  of  a  Reactive  Solid  by 
COMBUS  SCI  AND  TECH,  Vol.  2  (1970).  pp.  11-20. 
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